library(fixest)
library(modelsummary)
library(spmoran)
library(REndo)
library(WeightIt)
library(DirectEffects)
library(vtable)
library(spdep)
library(panelView)
library(ggplot2)
library(ggpubr)
library(viridis)
library(patchwork)
library(cragg)
library(mediation)
library(sensemakr)
library(robomit)
library(spatialreg)
library(rdrobust)

data <- readRDS("data_final_merged.rds")  %>% dplyr::select(NUTS_ID, Judet, Year, Month, p_ottoman, share, count, ps_1699, lon, lat, tavg, lograin, logelev, logslope, wheat_rainfed, log_large_rivers_distance,log_sea_lines_distance,route_density_2, vienna_distance, istanbul_distance,geometry)

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

data <- data %>%
  mutate(unique_id = paste0(NUTS_ID, "_", Year, "_", Month))

unique_data <- data %>%
  dplyr::distinct(unique_id, .keep_all = TRUE)

coords <- cbind(unique_data$lat, unique_data$lon)
nb <- knn2nb(knearneigh(coords, k = 4))  

nb_filtered <- lapply(1:length(nb), function(i) {
  neighbors <- nb[[i]]
  valid_neighbors <- neighbors[unique_data$NUTS_ID[neighbors] != unique_data$NUTS_ID[i]]
  
  if (length(valid_neighbors) == 0) {
    valid_neighbors <- neighbors[1]  
  }
  return(valid_neighbors)
})

attr(nb_filtered, "region.id") <- attr(nb, "region.id")
class(nb_filtered) <- "nb"

listw_filtered <- nb2listw(nb_filtered, style = "W")
listw<-listw_filtered 

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

model1 <- errorsarlm(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin +
                       log_large_rivers_distance + log_sea_lines_distance ,
                     data = data, listw = listw, method = "eigen")
model2 <- lagsarlm(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin +
                     log_large_rivers_distance + log_sea_lines_distance,
                   data = data, listw = listw, method = "eigen")
model3 <- lagsarlm(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin +
                     log_large_rivers_distance + log_sea_lines_distance,
                   data = data, listw = listw, Durbin = TRUE, method = "eigen")
model4 <- errorsarlm(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin +
                       log_large_rivers_distance + log_sea_lines_distance ,
                     data = data, listw = listw, Durbin = TRUE, method = "eigen")
model5<- sacsarlm(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin +
                    log_large_rivers_distance + log_sea_lines_distance,
                  data = data, listw = listw, method = "eigen")

msummary(list(model1, model2, model3,model4,model5), stars = c('.' = .1,'*' = .05, '**' = .01, '***' = .001), fmt = "%.3f")


